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Abstract 

This paper explores a possible technique for extending to multidimensional 
flows some of the upwind-differencing methods that have proved highly 
successful in the one-dimensional case. Attention here is concentrated on the 
two-dimensional case, and the flow domain is supposed to be divided into 
polygonal computational elements. Inside each element the flow is represented 
by a local superposition of elementary solutions consisting of plane waves not 
necessarily aligned with the element boundaries. 
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!• Introduction 


The recent survey of Woodward and Colella [16] shows that for one- 
dimensional gas dynamics there is an order of magnitude difference in 
effectiveness between sophisticated codes physically based on correct transfer 
of information, and simpler codes combining central differences with 
artificial viscosity. The sophisticated codes need much more computational 
work to update the solution at each mesh point, but this is far outweighed by 
their ability to capture discontinuities on a coarser mesh. For two- 
dimensional problems the difference in efficiency is far less marked, and for 
less violent flows than the ones they consider the advantage is likely to be 
reversed. 

The explanation is probably that the physics of one-dimensional flow is 
especially simple and well understood, and easy to imitate by numerical 
processes. Two-dimensional flows are more complex; in particular, acoustic 
waves can propagate in infinitely many directions rather than just two, and 
vorticity exists as a new phenomenon. Most extensions of upwind codes to two 
or more dimensions ignore these issues and advance the solution by 
splitting'’, that is to say, through a sequence of one-dimensional 
operators. For examples, see the survey by Woodward and Colella, also Sells 
[13] and Chakravarthy and Osher [1]. There are also what may be called "one 
and a half —dimensional" methods, in which the one— dimensional operators are 
interwoven, but the underlying physical model is still one involving wave 
propagation along the coordinate directions. This approach seems to yield 
some modest gains, as shown by Lytton [7] and Colella [2]. However, an 

observation is made in Section 4 which casts doubt on its real value. 
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If full advantage can be taken of upwinding techniques in two or more 
dimensions it is probably necessary to devise methods which take account of 
the actual directions in which information is propagated. The only results so 
far available for a method of this kind are those of Davis [4] . He assumes 
that the flow is locally dominated by a single shock wave whose unknown 
orientation may be deduced from the velocity field, or, in a later version of 
the code, from the pressure field (Davis, private communication). His method 
works very well on test problems where the flow is divided by shock waves into 
piecewise uniform regions. This is encouraging because it shows that a well 
chosen model of the flow can be used to numerical advantage. 

It has been conjectured that the way forward into two dimensions is 
blocked by the complexity of a 'two-dimensional Riemann solver', by which is 
meant an algorithm for computing the breakdown of initial conditions which are 
piecewise constant in two-dimensional cells. The solution of this problem 
close to the edge of a cell is straightforward, but secondary interactions 
near the corners are extremely difficult to compute. Even if a Riemann solver 
of this kind were computationally feasible, however, it would not be a 
satisfactory building block for two-dimensional calculations. It would, like 
the operator-splitting methods mentioned above, force the principal wave 
motions to take place normal to the cell boundaries. 

In the present work we avoid this difficulty by thinking of the data as 
piecewise linear rather than piecewise constant, and in Section 2 we interpret 
one-dimensional upwinding schemes in that light. The gradients in the data 
are used to construct a 'model flow' consisting of simple waves within each 
mesh interval. In Section 3 the corresponding simple wave solutions, 
propagating in arbitrary directions, are derived for the two-dimensional 
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equations. In Section 4 we propose model flows which can be fitted to any 
data which varies linearly in two dimensions, and in Section 5 we describe a 
strategy for constructing conservative differencing schemes by fitting such 
models to the data given at vertices of an irregular two-dimensional mesh. 
Section 6 contains observations on the possibility of extending the work to 
three-dimensional flow, and Section 7 comments on the type of advection scheme 
needed to complete the algorithm. 


2. Upwinding in One Dimension 

We begin by observing that one way to derive upwind schemes for the Euler 
equations in one dimension is to suppose that the flow in each mesh interval 
(i, i + 1) is a locally linear superposition of simple waves having the form 


w(x,t) = l r^x - \ t), 
k 


( 2 . 1 ) 


Here, w is the vector of unknowns, r^ is an eigenvector showing how the 


gradients due to the k c ^ wave are distributed over the components of w, 
is the amplitude of the k fc ^ wave, and X^ its speed. Any independent set 
of unknown variables w may be chosen, and the choice will not affect the 
values of a^, X^, but r^ will be different for each choice. The values of 


*1,2,3 are 


u - a, u, u + a 


( 2 . 2 ) 


1,2,3 


and the values of a 


are 
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— [ & p “ PaAvi] , Ap - — j Ap, — ^ [Ap + paAu] (2.3) 

2a a 2a 

where A(») = (»)^ + ^ - (*)^ and any local average values of p> a , u are 
valid. To achieve a conservative algorithm, two conditions are necessary. 
The eignevectors r^ must show the effects of the waves upon the conserved 
variables (p, pu, pe). In these variables r ^ , x^, r^ are given by 



1 


1 


1 


u - a 
h - ua 

- 2 

CM 

II 

si 

u + a 
h + ua 


2 i 2 

where h = a /(y - 1) + V2 u is the specific enthalpy. Also, the average 
values of p, a, u must now be chosen so that 

K \ £k ■ “l < 2 - 5 > 

where is the vector of flux quantities. For details, see [9,12]. 

This may be thought of as constructing, within each interval (i, i + 1) 
a local model of the flow. The model consists of elementary solutions of the 
Euler equations, linearized about a particular local average state. The model 
matches the observed data with respect to the spatial derivatives (or to be 
precise, with respect to the mesh differences). The time evolution of the 
model flow is readily predicted, and provides the information which is used to 
advance the global solution through one time step. 
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3. Elementary Solutions in TWo Dimensions 

In this section, as a necessary preliminary to the construction of two- 
dimensional models, we investigate the elementary solutions from which they 
may be built. Consider the Euler equations in primitive variables 
w = (p, u, v, p). 


2 

p. + up + vp + pa (u + v ) 
t *x r y x y 


u + uu + vu + — p 
t x y p r x 


v + uv + vv + — p 
t x y p y 


p_ + up + vp + p(u + v ) 
t x y x y 


Corresponding to Eq. (2.1) there are solutions 



to (3.1) of the form 


(3.1) 


w = a(0)r(0)(x cos 0 + y sin 0 - X(0)t) 


(3.2) 


where 0 is an arbitrary angle, and r(0) is an eigenvector. For acoustic 
waves and primitive variables it can be shown that 


r(0) 


2 

pa 

a cos 0 
a sin 0 


P 


(3.3) 


and that the wave speed is 



X(0) = u cos 0 + v sin 0 + a 


(3.4) 


Figure 1 shows that this speed corresponds to a wave front tangential to 
the Mach cone. We repeat here an observation from Roe [10]. Consider two 
Cartesian points having the same value of y, in a flow given by (3.2)- 
(3.4). An operator-spilt method will attempt to explain the difference in 
states as due to waves passing in the x-direction; it will compute 



P 


Pa 2 


pa 2 


0 


pa 2 


u 


a cos 0 


a 


0 


-a 

A 

V 

cc 

a sin 0 

= 1/2 (1 + cos 0) 

0 

+ sin 0 

a 

+ V2 (1 - cos 0) 

0 


p 


P 


P 


0 


P 


(3.5) 


where the RHS shows the eigenvectors of two one-dimensional acoustic waves, 
and a slip line. These spurious waves may not even travel in the proper 
direction and their inclusion in a numerical method can hardly be realistic. 
This criticism applies even to the "unsplit" algorithms of Colella [2] and 
Lytton [7]. Our goal in the next section is to construct local models of the 
flow by superposing simple waves whose orientation is not assumed in advance. 

Such a model cannot, however, be constructed purely out of acoustic waves 
since these are irrotational and the data may not be. There are two other 
fundamental flows which can be incorporated neatly into the model. One is a 
shear flow, which again has the general form of (3.2) but with 


r(0) 


0 

sin 0 
cos 0 


(3.6) 


0 
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X(0) = u cos 0 + v sin 0. (3.7) 


Another is solid-body rotation, or vorticity 


u = - V 2 uy v - V 2 wx . (3.8) 


There still remains an effect which is missing from the model, for all the 
flows above are isentropic. An entropy wave, across which pressure and 
velocity do not change, reveals itself in the primitive variables as a change 
of density. The general form is again (3.2) with X(0) given by (3.7), but 
with 

0 

0 

r = (3.9) 

0 

P 

independent of 0 . 

Another interesting fundamental solution (not directly used below) is 
obtained by superposing acoustic waves of the same strength with all possible 
propagation directions, i.e., by integrating (3.2) with respect to 0 from 
0 to 2rc with a(0) = a^. The result is 
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This solution would appear in the data as a region of uniform (isotropic) 
velocity divergence. However, the same data could be explained equally well 
by the passage of four plane waves 

w(0) + w(tt/2) + w(tt) + w(3tt/2) (3.11) 

where w(0) is given by (3.2). For numerical purposes the discrete 
representation by four plane waves is more amenable than the representation by 
one circular wave, and this is how a uniformly diverging flow would be dealt 
with in the model we develop below. However, it may be worth noting that any 
three equal waves separated by angles of 2ir/3 would also produce (locally) 
the same effect. 


4. The Discrete Models 

It is not obvious how the model flows of Section 2 should be generalized 
from one dimension to two. The chief difficulty is that whereas in one 
dimension there are just three types of elementary wave, in two dimensions 
there are infinitely many if we count all the possible orientations as 
distinct. In one dimension there is only one model that can be constructed, 
and it has three parameters which are the unknown wave strengths. Matching 
the model to the spatial gradients of the three data quantities p, u, p 
gives three simple linear equations whose solution is (2.3). In two 
dimensions the data will allow us to estimate gradients in two directions of 
four quantities, yielding eight items of information. Whatever model we 
choose must have eight free parameters, some of which may be wave amplitudes. 
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and the remainder will be orientation angles. If all the orientations are 
supposed to be known (aligned, for example, with the grid directions) we will 
again find easily solved linear equations for the amplitudes. However, 
because of the observation made concerning Eq. (3.5), we reject this approach, 
and require that at least some of the orientations be left unspecified. 
However, the equations which must be solved for the parameters then become 
nonlinear. If the free parameters are not judiciously chosen, no closed form 
solution may be possible, or the solution may not always be real-valued, or 
the solution may be computationally expensive. In such cases, the model will 
be useless. 

Two models, however, have been found whose parameters are given by simple 
real-valued expressions for all data. Each has, as its representation of the 
acoustic disturbances, a set of four orthogonal waves (Figure 2). One of the 
four will have an orientation angle in the range [±ir/4] and we take this as 
reference. Its orientation is 0, and its amplitude . The strength of the 
wave which moves in the opposite direction will be c^, and the waves which 
travel at right angles to these two have strengths . To this model we 
add an entropy wave with strength 0 and inclination <j>, so that the model 
now contains seven unknown parameters. 

To close the model we must introduce a fundamental solution incorporating 
vorticity, and it is only in this respect that the two models differ. In 
Model A we introduce a uniform vorticity w, and in Model B we introduce a 


shear flow such that 
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u = u Q (l + k(v Q x — Uq y)) 


v = v Q (l + k(v Q x - u q y)). 


(4.1) 


This is a special case of (3.6), (3.7) with tan 8 = -(u/v). We will 
first show the algebra for Model A, which is slightly simpler. 

To tidy up the equations we write dimensionless derivatives 


2 2 
P = p /pa P = p /pa 
x x y y 


U = u /a 
x x 


V = v /a 

X X 


U - u /a 

y y 


V = v /a 

y y 


(4.2) 


R * - p x /p \ ■ V 0 - 


By equating these to the sum of contributions produced by each component of 
the model, we find 


P ** a. cos 0 + a„ cos 0 - a. sin 0 - a. sin 0 

x 1 2 3 4 


(4.3a) 


P = a, sin 0 + a„ sin 0 + a. cos 8 + a. cos 0 

y 1 2 3 4 


(4.3b) 


2 2 2 2 
U = a, cos 0 - a. cos 0 + a 0 sin 0 - a, sin 0 

x 1 2 3 4 


(4.3c) 
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Uy = sin 0 cos 0 - sin 0 cos 0 - sin 0 cos 0 

+ sin 0 cos 0 - V 2 w/a (4.3d) 


V x = sin 0 cos 0 - sin 0 cos 0 - sin 0 cos 0 

+ sin 0 cos 0 + V; 2 w/a (4.3e) 

2 2 2 2 
Vy = a 1 sin 8 - sin 0 + a 3 cos 0 - cos 0 (4.3f) 

R x = otj cos 0 + a 2 cos 0 - sin 0 - sin 0 + 0 cos <|> (4.3g) 

Ry = a 1 sin 0 + sin 0 + a 3 cos 0 + cos 0 + 0 sin <{>. (4.3h) 

In these equations, the convention which distinguishes the contributions of 
otp is that the same angle 0 is used, but the sign of a is reversed. 

The eight equations can be solved quite easily. From (4.3d) and (4.3e) we 
obtain at once 

a) = a(V - U ) 
x y 

= (v - u ). (4.4) 

x y 
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whence 3, <{> . Next add (4.3d) and (4.3e) to obtain 


U + V = 2 sin 0 cos 0(a. 
y x 1 


a 2 - a 3 + a 4 ) 


(4.6) 


and subtract (4.3f) from (4.3c) 


U 

x 



(cos 2 0 - sin 2 0)(a 


a 2 - a 3 + a 4 ). 


(4.7) 


Dividing (4.6) by (4.7) yields 


tan 29 = 


U + V 

y x 
u 


X y 



(4.8) 


Since we have defined 1 9 j < tt/ 4 this result defines a unique orientation 
which is always real, coinciding, in fact, with the principal axis of the 
strain tensor. With 0 known, the remaining equations are linear. We write 
(4.3c) and (4.3f) as 


2 2 

U x = (c^ - a 2 )cos 0 + ( a 3 “ <*^)sin 0 

2 2 

Vy = (a^ - a 2 )sin 0 + (a^ - a^)cos 0 


and combine them to give 


U cos 2 0 - V sin 2 0 
x y 

°1 °2 


2 2 
cos 0 - sin 0 


(4.9) 
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This expression must be rewritten to avoid the possible singularity. Noting 
that 

9 9 

cos 0 =V2(1 + cos 20), sin 0 =1/2(1 “ cos 20) 

and that 

cos 20 = D/R (4.10) 

where 

R 2 = N 2 + D 2 (4.11) 

we find 

«1 - « 2 = l7 2 (U x + V y + R) (4.12) 

which is clearly always finite. By the same process we find 

a 3 - « 4 =V2 (U x + V y - R). (4.13) 

It can be shown that these expressions (4.12), (4.13) are proportional to the 
greatest and least straining rates experienced by the fluid. In these 
results, R must have the same sign as D, since 1 6 1 < tt/4 and so the RHS of 
Eq. (4.10) must be positive. For locally one-dimensional flow in the 
x(resp y) direction, R will equal U x (resp -V y ) and V y (resp U x ) will be 
zero. Eqs. (4.12), (4.13) will give the correct one-dimensional results. 

That is, “ a 2 = U x (resp 0), and = 0(resp V^,). It is interesting 

that 

<*1 " « 2 + a 3 " a 4 = U x + V (4 ' 14) 

The LHS is the total strength of the acoustic waves (the minus signs appear 
because of our conventions about a and 0) and the RHS is the velocity 
divergence. Compare the result in Eq. (3.10). 
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The last step is to combine (4.3a), (4.3b) to give 

a. + a„ = P cos G + P sin 9 (4.15) 

1 2 x y 

a 0 + a, = P cos 9 - P sin 0 (4.16) 

3 4 y x 

and then the a's follow from (4.12), (4.13). 

A. remarkable identity concerning the wave strengths is the following; 

a l + °2 + a 3 + a 4 + ^ 

2 

“ V2 (<*j + a 2 ^ + ^2 “ a 2 + V2 ( a 3 + a 4)^ +V2( a 3 " <*4) + — *r 

4a 

= 1/2 (P cos 0 + P sin 0) 2 + i (U + V + R) 2 
x y 8 x y 

2 

+ 1/2 (P cos 0 - P sin 9) 2 + i (U + V - R) 2 + -Ar 
L y x 8 x y 4a 2 

2 

= V 2 (pJ + P?) + V4 (U + V ) 2 + 1/4 R 2 + ~ ~2 

y 3 4a 

- V 2 (p£ + Py) + V4 (U x + v y ) 2 +1/4 (U x - v y ) 2 +1/4 (u y + V x ) 2 +1/4 (V x - U y ) 2 

= 1/2 [p 2 + Py + Ux + U y + V X + V yl* (4 * 17) 

Both ends of this chain are expressions representing some overall strength 
of the disturbance (excluding entropy effects which add another simple term). 
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The analysis of Model B is almost identical. The equations are altered by 
replacing the vorticity terms with the shear terms from (4.1) thus 


P = a- cos 0 + a 0 cos 0 - a 0 sin 0 - a, sin 0 

x 1 2 3 4 


(4.18a) 


P = a, sin 0 + a 0 sin 0 + a„ cos 0 + a, cos 0 

y 1 2 3 4 


(4.18b) 


U x = cos ^ 0 - cos^ 0 + sin^ 0 - sin^ 0 — ku^ Vq (4.18c) 


U = a. sin 0 cos 0 - sin 0 cos 0 - cu sin 0 cos 0 + a, sin 0 cos 0 - ku^ 
y I 2 3 4 0 

(4 . 18d) 
2 

V = a sin 0 cos 0 - a n sin 0 cos 0 - a 0 sin 0 cos 0 + a. sin 0 cos 0 + kv_ 
xl 2 3 4 0 


(4 . 18e) 


V y - sin^ 0 - sin^ 0 + ot^ cos^ 0 - cos^ 0 - kuQ Vq (4.18f) 


R = a. cos 0 + a 0 cos 0 - a 0 sin 0 - a, sin 0+8 cos 6 

x 1 2 3 4 


(4 . 18g) 


R y = sin 0 + sin 0 + cos 0 + cos 0+8 sin <|>. 


(4 . 18h) 


The solution for 8, <f> is identical. Equations (4.18d), (4.18e) give 


V - U 

k = js l 

2 2 

u 0 + v 0 


(4.19) 


The expression for 0 in this case is 
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tan 20 


» T + V x * k ( u 0 ' v 0 ) 
u * - V y - 2ku 0 v 0 


(4.20) 


which can be rewritten, using (4.19), and setting v/u - tan 6, as 


tan 29 = 


u y + V x + <V X - U y )cos 26 _ N 
= u - (V - U )sin 26 D 


(4.21) 


Note that for irrotational flow, (4.21) agrees with (4.8). Again we 
introduce R, such that R 2 = N 2 + D 2 , and having the same sign as D. In 
terms of this new R, we still have 


= v 2 (U x + v y + R) 

(4.22) 

= v 2 (U x + v y - R) 

(4.23) 


and Eqs. (4.15), (4.16) are unaffected. 

For Model B there seems to be no simple analogue of Eq. (4.17). 

Otherwise, the difference between the two models is that Model B involves 

computing slightly more expensive expressions for N and D, but may be able 
t 0 fit itself to a greater variety of flows. Both models have the property 
that if the data is locally one-dimensional in any direction then waves will 
be predicted which are exactly those predicted by a one-dimensional linear 

Riemann solver aligned in that direction (rather than with the coordinate 

axis). However, Model B can simultaneously recognize a shear flow in some 
other direction. Neither model, however, could correctly recognize both shock 
waves of a colliding pair, unless these happened to be perpendicular. It 
would appear that any model flow must be a compromise between simplicity and 



-17- 


generality. Simple models will generally be invalid at isolated points, and 
reliance must then be placed on conservation. It is to this aspect that the 
next section is devoted. 


5 . Conservation Properties 

To create an algorithm capable of capturing shock waves, we must ensure 
that it is conservative. For present purposes, the most convenient definition 
of a conservative algorithm is that when it operates for one time step, the 
conserved quantities (mass, momentum, and energy) present within the 
computational domain are only changed because of events occurring on the 
boundaries of the domain. We will first set out a strategy which guarantees 
this. Then we will relate the results of previous sections to that strategy. 

Suppose that the computational domain is tessellated into arbitrary 
polygons (see Figure 3). Usually these would be quadrilaterals or triangles, 
and the formulae given below will then be very simple. However, we treat the 
general case to show that exceptional meshes create no difficulty, at least 
with regard to conservation. Consider, then, an arbitrary cell with vertices 
Vj, ^2***^n anc * note that the area of the cell may be written 

4A ' E x ( £i + i - £t-i> <5.0 

where _r^ is the position vector of the i fc ^ vertex, and the counting is 
cyclic and anticlockwise. Eq. (5.1) is proved by observing that the terms in 
the summation occur in equal pairs, and that every term x r^ + ^ is twice 
the area of a triangle V i 0V i+1 where 0 is an arbitrary origin. 

Rearrangement of terms in (5.1) leads to two alternative expressions 
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2a - l x i<y 1+ i - y ± -i> 


“ y i (x i+i " x i-i ) * 


(5.2) 


(5.3) 


Simple alterations of these formulae allow us to estimate the gradients 
within a cell of any quantity q which is defined at the vertices. Thus 


2A - I 1i<y 1+1 - Vt-P 

2A If ‘ -l ‘ Vi' 


and it can be seen that these estimates are exact whenever q is a linear 
function (q = mx + ny) . 

Now suppose that the quantities stored at the vertices are the variables 
defining flow of an ideal gas according to the Euler equations, written in 
conservation form as 


w + F + G = 0. 
— t —x — y 


(5.6) 


Then an estimate for w. , averaged over the cell, is 

— z 

“Ee * -I - G i <x i+ i - 


(5.7) 


An alternative way to obtain this formula is to integrate the passage of flux 
across the cell boundary, using the trapezium rule. We have followed this 
present derivation because the formulas (5.4), (5.5) are also useful for 

estimating the gradients from which, in Section 4, the local flow model was 


deduced. 
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The quantity , multiplied by a finite time step At represents the 
local accumulation of the conserved quantities. The solution can be advanced 
one time step by adding this change to the quantities stored at the 
vertices. The increments may be distributed equally or unequally to the 

vertices concerned. An equal distribution would, if applied on a regular 
rectangular mesh, reduce to a central-differencing scheme of the kind that can 
be allied with Runge-Kutta schemes [5]. An unequal distribution of 

increments, where the weights are obtained from the Jacobian matrices 3F/3w 
and 3G/3w, has been used by Ni [8] to obtain an integration which is 
equivalent to Lax-Wendrof f . The present work is intended for use with a 
scheme in which the increments are distributed with more regard to the 
'upwind' direction of each wave. Meanwhile, we prove that any distribution 
will lead to a conservative algorithm. 

The total change of conserved quantities, within the computational domain, 
is obtained by summing (5.7) over all cells. A typical vertex Vj in the 
interior of the domain, contributes to this sum through all the cells which 
meet there. Its total contribution is, in fact 

% = ["Ij I A y + I Ax] (5.8) 

where the Ax, Ay are the adjacent chords of each cell meeting at Vj (see 
Figure 3). But since the union of these chords is a closed polygon Aw. = 0. 
Since this argument applies equally to all interior vertices, the sum of 
conserved quantities changes only due to events on the boundary, and this is 
what we require. 
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Next we demonstrate how the estimated total increment (5«7) may be 
decomposed into contributions due to each wave system. It has not been found 
possible to do this by any direct extension of the analysis in Section 4. 
When the spatial changes are large, there seems to be no simple choice of mean 
values which allows a tidy analysis of the flux gradients. Instead, we 
directly analyze the temporal changes inside each cell to produce a 
decomposition which is conservative but not unique. Uniqueness is imposed by 
incorporating results from Section 4. 

First, observe that the time derivative of w due to the passage of a 
plane wave is the product of the amplitude and wave speed multiplied by an 
eigenvector which describes the effect of that wave on the conserved 
variables. For an acoustic wave inclined at an angle 0, that eigenvector is 


r 

—a 


P 

pu + pa cos 0 
pv + pa sin 0 

ph + pa(u cos 0 + v sin 0) 


(5.9) 


and for an entropy wave at any angle it is 


r “ 


P 

pu 

pv 

V2 p(u 2 + V 2 ) 


(5.10) 


For a shear wave it is 
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0 

-pa sin 9 
pa cos 0 

pa(v cos 0 - u sin 0) 


( 5 . 11 ) 


However, the shear wave included in Model B has zero speed (i.e., it is a 
steady solution of the Euler equations) so that the term involving r makes 
no contribution to w^_ . In this respect Model B is somewhat simpler than 
Model A, because the uniform vorticity does contribute to w t> in a way which 
is derived below. Introduce the notation 


Aj = = a j( u cos 8 + v sin 0 + a) (5.12a) 

A£ = a 2 ^2 = a 2^ u cos ® + v s * n ® " a ) ( 5 . 12 b) 

A3 = 03 X3 = o^C-u sin 0 + v cos 0 + a) (5.12c) 

A 4 = % X 4 " a 4 (- u sin 0 + v cos 0 - a) ( 5 . 12 d) 


A^ = = p(u cos <J> + v sin <)>). 


( 5 . 12e) 


Our strategy is to compute the {a^} within each cell in such a way that 
the total effect of all the disturbances in that cell will produce the correct 
conservative value of w t » First, though, it must be checked that the model 
does contain all the effects contributing to w fc . Therefore, we evaluate 


i =4 

Jj A i — ai + A 5 V 


( 5 . 13 ) 
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As an example, the terms contributing to are 


p[Aj + + A^] = pctj(u cos 0 + v sin 0a + a) 


+ pa^(u cos 0 + v sin 0 - a) + pa^(v cos 0 - u sin 0 + a) 


+ pa^(v cos 0 — u sin 0 - a) + pf3(u cos <|> + v sin <}>) 


= p(a^ + c^Ku cos 0 + v sin 0) + p(a^ - <x^)a 


+ pCa^ + a^)(v cos 0 - u sin 0) + P(Og - a ^) a 


+ p3(u cos <|) + v sin $) . 


Substituting the results of Section 4 into this expression, we obtain 


p[A^ + + A^ + A^ + A,.] = p[P x cos 0 + P y sin 0](u cos 0 + v sin 0) 


+ 



cos 0 - P x sin 0](v cos 0 - u sin 0) 


+ V 2 p[U x + V y + R]a +V2P[\ + V y - R]a 


+ pu[R x - P x ] + pv[R y - P y ] = pa[U x + V y ] + puR x + pvR y 


or, in terms of the dimensional gradients 



-23- 


p t A 1 + A 2 + A 3 + A 4 + a 5 J ~ pf u x + v y ] + up x + vp y = Pj; . (5.14) 

This calculation, which is valid for Models A or B, checks the algebra and 

confirms the completeness of the model. Checking the other components of w 
is tedious, but necessary. It reveals that Model B supplies all the terms of 
— t from the expression (5.13), but that when this expression is used to 
calculate the effects of the acoustic and entropy waves in Model A, there is a 
surplus in the expression for (pu) fc amounting to 1/2 pvw, and a shortfall in 

the expression for (pv) t of V2 Puw. These terms represent the effects of 

convected vorticity. The expression for (pe) t turns out to be correct. 

Therefore, we write 


p t = P(A 1 + A 2 + A 3 + A 4 + V 


(5.15a) 


(pu) t = pAj (u + a cos 0) + pA 2 (u - a cos 0) + pA 3 (u - a sin 0) 


+ pA^(u + a sin 0) + pA $ u - V2 P™ (5.15b) 


(pv) t = pA x (v + a sin 0) + pA 2 (v - a sin 0) + pA 3 (v + a cos 0) 


+ pA 4 (v - a cos 0) + pA 5 v + l / 2 P uaj (5.15c) 


(pe) t = pAj(h + au cos 0 + av sin 0) + pA 2 (h - au cos 0 - av sin 0) 


+ pA 3 (h - au sin 0 + av cos 0) + pA 4 (h + au sin 0 - av cos 0) 


+ 1/2 P a 5 (u 2 + v 2 ). 


(5.1 5d ) 
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where w = 0 for Model B. 

These equations are the two-dimensional analogue of Eqs. (2.5)* In each 
case we try to ensure that the changes of conserved variables predicted by the 
model are correct. Here, we assume that the LHS of each equation is obtained 
from the conservative formula (5.7) for each cell. Then we treat (5.15) as a 
set of conditions to be identically (not just approximately) satisfied by the 
{A ± } and by 0, 4>, (d. Since there are only four conditions for eight 
unknowns, the remaining information must be supplied from elsewhere. It seems 
natural to take the values of 0, <(>, w from Section 4. Conditions (5.15) are 
then an incomplete set of linear equations for the {A^}, which may be 
partially analyzed as follows. We obtain at once 

(pu) ^_ - up^ + pvw = pa cos 0(A^ - A^) - pa sin 0(A^ - A^) (5.16a) 

(pv) ^ - vp - puw = pa sin 0(A^ - A^) + pa cos 0(A^ “ A^) (5.16b) 

and hence (A^ - k^) , ( A - A^) . Substituting these results into (5.15d) 

yields 

2 

t A 5 = (h - u 2 - v 2 )p t + u(pu) t + v(pv) t - (pe) t . (5.17) 

If the changes are small, so that (•) may be treated as a derivative, 

these equations simplify considerably, offering more insight into the models. 

(A^ - A^) = u^ cos 0 + v fc sin 0 + V 2 (v cos ® " u sin 0)w (5.18) 


(A^ - A^) = v t cos 0 - u fc sin 0 - V 2 ( u cos 0 + v sin 0)w (5.19) 
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a 5 - (a 2 P t - P t )/pa 2 (5.20) 

A 1 + A 2 + A 3 + \ = P t « (5.21) 

Non-conservative schemes could use these simpler conditions; a fully 

conservative scheme would have to satisfy (5.15). 

One more condition is needed on the {a^. In view of the symmetry of the 

results so far, we would like an expression for A, + A„ - A - A (which 

12 3 4 

need not derive from the conservation form). By changing some signs in the 
analysis leading to (5.14) we find 

A ^ + A 2 ~ A g “ = p x ( u cos 20 + v sin 20) 

+ P y (u sin 20 - v cos 20) + aR. (5.22) 

It may be shown that the RHS does not, in general, vanish when the data 

are taken from a steady flow. One might suppose that it should, since then 

the {A i } would all be zero, and either the strength or the speed of every 

wave would be zero. Instead of this, the models represent steady flow by a 

state of equilibrium between finite waves, such that A *» A = -A = -A 

1 2 3 4 * 

and A,. = 0. 

We have now generated a conservative model of the flow, in which the 
effects of the various components are given by (5.15). The parameters of this 
decomposition ({a ± }, 0, $, to) are found from (4.4), (4.5), (4.8), or (4.21), 
(5.16), (5.17), and (5.22). Any consistent choice of local average values for 
P, u, v, a, h in these equations will be valid, and will not affect the 
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conservation property. It may be asked, though, whether there are particular 
average values, similar to those which appear in the one-dimensional theory 
[9,12], bestowing special 'shock-recognition' properties. However, this 
question raises unsolved problems about the sort of 'captured shock structure' 
that is possible in two-dimensional flow, and will not be discussed here. 


6. Extension to Three Dimensions 

No detailed formulae have been worked out for the three-dimensional 
case. However, merely counting the degrees of freedom makes it plausible that 
analogous models could be constructed. Data for the three-dimensional 
unsteady Euler equations would consist of five variables, so there would be 
fifteen gradients to be accounted for by the model. If the acoustic 
disturbances are again to be represented by a set of orthogonal plane waves 
(like an expanding cube) there will be six wave amplitudes and three angles 
involved (two angles to orient one wave, one angle to orient its neighbors). 
An entropy wave with one amplitude and two angles will bring the number of 
parameters up to twelve. The remaining three are available to represent 
rotational effects. The analogue of Model A would contain three independent 
vorticity components. The analogue of Model B could contain a shear flow 

q = qjl + k^VQ x - u Q y) + k 2 (w Q y - v Q z) + k 3 (u Q z - w Q x) ] (6.1) 

which is, like (5.1), a steady solution to the Euler equations. However, the 
three shear components which it contains are not all independent, since all 
take place in a parallel flow, and one of the can be dropped with no loss 
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of generality. To complete the set of fifteen parameters one might add the 
flow 

u = u Q = 0)(w 0 y - v Q z) 

v = v 0 + »(u 0 z - w Q x) (6.2) 

w = w Q + io(Vq x - u Q y). 

This is also a steady solution of the Euler equations and represents a 
swirling flow in which the vorticity is parallel to the streamline (uq> vq, 
wq ) • Again there is a computational advantage to Model B in that some of its 
components are steady flows whose contribution to the time-marching process 
are identically zero. In fact, an analogue of Model B can be worked out for 

any number of space dimensions d, and the description of arbitrary data is 

reduced to the description of (2d + 1) non-linear scalar advection problems. 

There are, however, geometrical difficulties which appear in three 
dimensions, if the partition of space is made into volumes whose facets have 
more than three sides. Most finite-volume schemes employ computational cells 
which are hexahedral, with quadrilateral faces specified by four vertices not 
normally lying in one plane. The boundary surfaces of the cells must be 
unambiguously defined so that the cells fill the computational space without 
overlaps or voids. This can be done by folding each surface into two 
triangles [6], or by choosing a particular doubly-ruled surface to cover each 
face [3]. Once this has been done, any consistent formula for the volume can 
also be used to estimate the flux divergence, as in (5.7), but the consistent 
formulae are disconcertingly lengthy [6,3]. The effect of simpler formulae on 
the computational accuracy awaits investigation. 
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7 . Advection Schemes 

The eventual goal of this work is to create an algorithm for 
multidimensional gas dynamics which will enjoy the same degree of success 
already obtained by upwind schemes in one dimension. In this paper we have 

addressed (at most) half the problem, showing how arbitrary disturbances in 

the data can be replaced by locally equivalent sets of plane waves and/ or 
vorticity. To march forward in time, we need to apply to each wave some 
numerical advection scheme. We may hope that such schemes can be based on 

schemes for scalar problems, as has happened in the one— dimensional case. 

Also, we may anticipate that such schemes will show many of the typical 
features of successful one-dimensional schemes, such as asymmetric support and 
non-linear limiters [15]. However, the theory even of scalar advection 
algorithms in many dimensions is only in its infancy. Roe and Baines [11] 
present a criterion designed to avoid overshoots and describe a scheme which 
meets it. Smolarkiewicz [14] describes another distinctive, but related, 
approach. The next (rather large) step in the investigation reported here 
will be to experiment with these and other algorithms in the present context. 


8 . Conclusions 

We have pointed out that the extension of upwind differencing schemes to 
more than one space dimension cannot be accomplished by operator splitting 
methods without losing the desirable property of recognizing data due to a 
simple wave. To construct "genuinely two-dimensional" schemes we propose 

model flows, composed of elementary solutions to the two-dimensional 
equations. These model flows are such that they can be matched to arbitrary 
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data which varies linearly in some small region. The acoustic part of the 
flow is modelled by four orthogonal plane waves whose orientation is matched 
to the gradients in the data. Variation of entropy is represented by a single 
plane wave, and rotational effects either by uniform vorticy or by a parallel 
shearing motion. We show that the parameters of the model can be evaluated in 
such a way that a time- marching algorithm can be made exactly conservative. 
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Figure 1 . An acoustic wave front. 
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